## "Public Preferences for International Law Compliance: 
##  Respecting Legal Obligations or Conforming to Common Practices?"
##
##
## Reproduce: 
##      - Figure 1 
##
## Input:  
##      - data/data_cleaned.csv
##
## Output: 
##      - f1.pdf
##
## Saki Kuzushima
## 01/06/2023

source("code/functions_f1.R")
data <- read.csv(file="data/data_cleaned.csv")


# output table
item <- c("surname3", 
          "whale1",
          "hate1",
          "death1")
treat_group <- c("intl", "const", "oecd")
res <- expand.grid(item=item, treat=treat_group, stringsAsFactors=F)
res$lower <- 0
res$point <- 0
res$upper <- 0
res$p <- 0
res$se <- 0
res_simple <- res

# block rand. estimator
for (i in seq(1, nrow(res))){
  treatment <- paste0("treat_", res[i,"treat"])
  item <- res[i, "item"]
  res[i, c("lower", "point", "upper", "p", "se")] <- 
    diffmean_block(data=data, treatment=treatment, item=item)
}


# multiple testing correction
# BH
res$p_bh <- p.adjust(res$p, method="BH")

# ash
ash_fit <- ash(betahat=res$point, 
               sebetahat=res$se, 
               prior="uniform", 
               mixcompdist="unif",
               grange=c(-5, 5),
               outputlevel=2)
ash_result <- ash_fit$result
ash_ci <- ashci(ash_fit)
res_ash <- data.frame(item=res$item,
                      treat=res$treat,
                      lower=ash_ci[,1], 
                      point=ash_result[,"PosteriorMean"],
                      upper=ash_ci[,2])

# fiugre 
x_index <- c(1.25, 1.75, 3, 4.25, 4.75, 6, 7.25, 7.75, 9, 10.25, 10.75) 
point_intl <- c(res$point[1], res_ash$point[1], NA, 
                res$point[2], res_ash$point[2], NA, 
                res$point[3], res_ash$point[3], NA, 
                res$point[4], res_ash$point[4])
upper_intl <- c(res$upper[1], res_ash$upper[1], NA, 
                res$upper[2], res_ash$upper[2], NA, 
                res$upper[3], res_ash$upper[3], NA, 
                res$upper[4], res_ash$upper[4])
lower_intl <- c(res$lower[1], res_ash$lower[1], NA, 
                res$lower[2], res_ash$lower[2], NA, 
                res$lower[3], res_ash$lower[3], NA, 
                res$lower[4], res_ash$lower[4])


point_const <- c(res$point[5], res_ash$point[5], NA, 
                 res$point[6], res_ash$point[6], NA, 
                 res$point[7], res_ash$point[7], NA, 
                 res$point[8], res_ash$point[8])
upper_const <- c(res$upper[5], res_ash$upper[5], NA, 
                 res$upper[6], res_ash$upper[6], NA, 
                 res$upper[7], res_ash$upper[7], NA, 
                 res$upper[8], res_ash$upper[8])
lower_const <- c(res$lower[5], res_ash$lower[5], NA, 
                 res$lower[6], res_ash$lower[6], NA, 
                 res$lower[7], res_ash$lower[7], NA, 
                 res$lower[8], res_ash$lower[8])

point_oecd <- c(res$point[9], res_ash$point[9], NA, 
                res$point[10], res_ash$point[10], NA, 
                res$point[11], res_ash$point[11], NA, 
                res$point[12], res_ash$point[12])

upper_oecd <- c(res$upper[9], res_ash$upper[9], NA, 
                res$upper[10], res_ash$upper[10], NA, 
                res$upper[11], res_ash$upper[11], NA, 
                res$upper[12], res_ash$upper[12])

lower_oecd <- c(res$lower[9], res_ash$lower[9], NA, 
                res$lower[10], res_ash$lower[10], NA, 
                res$lower[11], res_ash$lower[11], NA, 
                res$lower[12], res_ash$lower[12])

# Style memo:
# Point style: hollow for uncorrected, filled for corrected ones
# Note that pch=18 (filled diamond) looks smaller than others. 
# So the cex for pch=18 is set to 2.0; others 1.5
# CI color: grey for non-significant ones; Black for significance
point_style <- c(1,16,NA,2,17, NA,0,15, NA,5,18)
pdf("results/main/f1.pdf")
par(mfrow=c(3,1), mar=c(1,4,1,1)+0.1, oma=c(2,1,3,1)+0.1)
# Row1: International Law
plot(x=x_index, y=point_intl, type="n",
     xlim=c(0.1, 11.9), ylim=c(-0.5, 0.5),
     xaxs="i", yaxs="i",
     yaxt="n", xaxt="n",
     xlab="", ylab="International Law", cex.lab=1.25)
abline(h=0, col="red")
abline(v=3, lty=3)
abline(v=6, lty=3)
abline(v=9, lty=3)
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
axis(side=3, at=c(1.5, 4.5, 7.5, 10.5), labels=c("Surname", "Whale", "Hate Speech", "Death Penalty"),
     tick=F,  line=0, cex.axis=1.25)
points(x=x_index[c(1,4,7,10)], 
       y=point_intl[c(1,4,7,10)],
       pch=point_style[c(1,4,7,10)],
       cex=1.5,
       col="grey50") 
points(x=x_index[c(2,5,8,11)], 
       y=point_intl[c(2,5,8,11)],
       pch=point_style[c(2,5,8,11)],
       cex=c(1.5, 1.5, 1.5, 2),
       col="grey50") 
points(x=x_index[7], 
       y=point_intl[7],
       pch=point_style[7],
       cex=1.5,
       col="black") 
segments(x0=x_index[c(1,4,7,10)], 
         x1=x_index[c(1,4,7,10)],
         y0=lower_intl[c(1,4,7,10)],
         y1=upper_intl[c(1,4,7,10)], col="grey50") 
segments(x0=x_index[c(2,5,8,11)], 
         x1=x_index[c(2,5,8,11)],
         y0=lower_intl[c(2,5,8,11)],
         y1=upper_intl[c(2,5,8,11)], col="grey50") 
segments(x0=x_index[7], 
         x1=x_index[7],
         y0=lower_intl[7],
         y1=upper_intl[7], col="black") 
# Row 2: constitution
plot(x=x_index, y=point_const, type="n",
     xlim=c(0.1, 11.9), ylim=c(-0.5, 0.5),
     xaxs="i", yaxs="i",
     yaxt="n", xaxt="n",
     xlab="", ylab="Constitution", cex.lab=1.25)
abline(h=0, col="red")
abline(v=3, lty=3)
abline(v=6, lty=3)
abline(v=9, lty=3)
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
points(x=x_index[c(1,4,7,10)], 
       y=point_const[c(1,4,7,10)],
       pch=point_style[c(1,4,7,10)],
       cex=1.5,
       col="grey50") 
points(x=x_index[c(2,5,8,11)], 
       y=point_const[c(2,5,8,11)],
       pch=point_style[c(2,5,8,11)],
       cex=c(1.5, 1.5, 1.5, 2),
       col="grey50") 
points(x=x_index[c(1,4,5)], 
       y=point_const[c(1,4,5)],
       pch=point_style[c(1,4,5)],
       cex=1.5,
       col="black") 
segments(x0=x_index[c(1,4,7,10)], 
         x1=x_index[c(1,4,7,10)],
         y0=lower_const[c(1,4,7,10)],
         y1=upper_const[c(1,4,7,10)], col="grey50") 
segments(x0=x_index[c(2,5,8,11)], 
         x1=x_index[c(2,5,8,11)],
         y0=lower_const[c(2,5,8,11)],
         y1=upper_const[c(2,5,8,11)], col="grey50") 
segments(x0=x_index[c(1,4,5)], 
         x1=x_index[c(1,4,5)],
         y0=lower_const[c(1,4,5)],
         y1=upper_const[c(1,4,5)], col="black") 
text(x=5, y=0.25, labels="BH"~symbol("\326"))
# Row 3: Common practices 
plot(x=x_index, y=point_oecd, type="n",
     xlim=c(0.1, 11.9), ylim=c(-0.5, 0.5),
     xaxs="i", yaxs="i",
     yaxt="n", xaxt="n",
     xlab="", ylab="Practices", cex.lab=1.25)
abline(h=0, col="red")
abline(v=3, lty=3)
abline(v=6, lty=3)
abline(v=9, lty=3)
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
points(x=x_index[c(1,4,7,10)], 
       y=point_oecd[c(1,4,7,10)],
       pch=point_style[c(1,4,7,10)],
       cex=1.5,
       col="grey50") 
points(x=x_index[c(2,5,8,11)], 
       y=point_oecd[c(2,5,8,11)],
       pch=point_style[c(2,5,8,11)],
       cex=c(1.5, 1.5, 1.5, 2),
       col="grey50") 
segments(x0=x_index[c(1,4,7,10)], 
         x1=x_index[c(1,4,7,10)],
         y0=lower_oecd[c(1,4,7,10)],
         y1=upper_oecd[c(1,4,7,10)], col="grey50") 
segments(x0=x_index[c(2,5,8,11)], 
         x1=x_index[c(2,5,8,11)],
         y0=lower_oecd[c(2,5,8,11)],
         y1=upper_oecd[c(2,5,8,11)], col="grey50") 
dev.off()
